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ABSTRACT 


Control of a modem submarine is a multi-dimensional problem coupling initial 
stability, hydrodynamic and control system response. The loss of stability at moderate 
to high speeds is examined using a nonlinear Hopf bifurcation analysis. Complete linear 
state feedback is used for demonstration purposes for depth control at level attitude and 
for a fixed nominal speed. Control time constant, nominal and actual speeds, metacentric 
height, and stem to bow plane ratio are used as the main bifurcation parameters. A 
complete local bifurcation mapping provides a systematic method for evaluating the 
bounds of controllability for control system design parameters for a submarine with a 
given set of hydrodynamic coefficients. The submarine and its potential design 
modifications are then verified with a nonlinear simulation program. 
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closed loop dynamics matrix for the linearized 
system 

dummy independent variable 

bow plane to stern plane deflection ratio 

(or, control surface coordination gain) 

local beam of the hull 

vehicle buoyancy 

quadratic drag coefficient 

coupled heave and pitch terms 

cross flow drag terms 

bow plane deflection 

stern plane deflection 

saturation value of 5 

vehicle mass moment of inertia 

controller gains in fi,w,q, and z, respectively 

cubic stability coefficient 

vehicle mass 

pitch moment 

derivative of M with respect to a 

pitch rate 

vehicle pitch angle 

polar coordinates of transformed reduced system 

matrix of eigenvectors of A 

time constant 

vehicle forward speed 

nominal forward speed 

heave velocity 

vehicle weight 

state variable vector 

body fixed coordinates of vehicle center of buoyancy 
body fixed coordinates of vehicle center of gravity 
state variables vector in canonical form 
deviation off the nominal depth 
vehicle metacentric height 
heave force 

derivative of Z with respect to a 
frequency at the bifurcation point 
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I. INTRODUCTION 


The fundamental goal of submarine control is to reach and 
maintain ordered depth. Any design that does not meet this 
goal, in the face of the inherent complexities, is not overly 
useful as a practical vessel. Current evaluation schemes 
involve extensive model testing such as those done on the 
DARPA SUBOFF Model (DTRC Model 5470) [ref. 8]. This is an 
expensive and time consuming evaluation method. The goal is 
to develop an analytic method to determine the stability of a 
potential design. Much work has been done on depth control 
and modelling of submarines in the vertical plane [ref's. 1, 
2, & 3] but this work has asstuned the existence of a stable 
platform. 

The stability of a design will have a significant impact 
on its responsiveness. A vehicle with a large margin to 
instability will not be very responsive. The problem becomes 
one of determining how close to the margins we can get without 
a total loss of stability. Nonlinear dynamics and chaos 
provides us with the tools for solving this problem [ref's 9, 
& 10] . By expanding the scope of our research to include 
nonlinear terms we are able to define the limits of stability 
and therefore the margins. At the Naval Postgraduate School 
there has been extensive work on defining the nature of the 
instabilities that are encountered in the control of 
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submersibles in path keeping and the dive plane [ref's 4, 5, 
& 6 ] . 

At low forward velocities a sxibmersible using stern planes 
for attitude and depth control can experience a loss of 
stability in the form of stem planes reversal. This is a 
pitchfork bifurcation that can be predicted [ref. 7] and can 
be accounted for in the control design. The purpose of this 
thesis is to develop a program for finding the limits of 
stability for a submarine at moderate and high speeds. Once 
these limits are mapped then the nature of the loss of 
stability must then be determined. For this we use a Hopf 
bifurcation analysis along with a nonlinear simulation 
program. Finally, after the limits are determined we are able 
to define the control system design parameters and evaluate 
the controllability of the design. 
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II. PR0BL8M POSMULATION 


A. SQUATIONS OF NOTIOK 

By restricting the submersible's motion to the vertical, 
or dive plane, the motion can be modeled by the coupled 
nonlinear equations for pitch and heave. With a body fixed 
coordinate frame at the vehicles geometric center, vre can 
express Newton's equations of motion as 


wiw - Uq - 2(;q 


2 _ 


x^) = Z^g + Z^w + ZgUq + Z„Uw 
TAIL |v-xq| 

(W-B)cos0 + U^iZ^bg + Z^b^,) , 


( 2 . 1 ) 


and 


Iy<^ + mz^wq - mx^iw - Uq) = M^q + + M^Uq + M^,Uw 

/ NOSE , , lw~ XCf) ^ , 

b{x) - iS^xdx 

TAIL \W - Xa 


TAIL |w - xq| (2.2) 

- {XgW - XgB) COS 0 - {ZqW - ZgB) sin 0 


iOr simplicity we will assume that the drag coefficient, C^, 
is constant over the length of the submersible. This does not 
have a significant effect on the results. 
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The pitch rate and depth rate for a submersible is given 
by 


± = -UsinQ + WCOS0 , (2.4) 

respectively. In equations (2.3) and (2.4) 6 is the pitch 
angle of the vehicle as measured from the horizontal plane. 
Figure (2.1) shows the geometry and definitions for most of 
the symbols. We assume that forward velocity, U, is 
maintained constant during maneuvers by an appropriate use of 
the propulsion system. 

B. CONTROL LAW 

Equations (2.1) through (2.4) can be written as a set of 
nonlinear differential equations in the form. 


w = a^^Uw + a^2^g + a^jZ^gSinS + 

+ d„{w.q) + Cl ( w, q) , 


( 2 . 6 ) 



g = a2^Uw + a^zUg + ajjZagSine + + b^^U^b^, 

+ dq{w,g) + C2(v, g) , 


(2.7) 


i = -C7sin0 + wcos6 , 


( 2 . 0 ) 


where, 

= im-Z^) (ly-M^) -{mx^+Z^) {mx^+M^) , 
a^J^^ = ily-M^) Z„+ imxa-^-Z^) M„ , 
ai2l>v = (ly-M^) (m+Zg)+{mXs+Z^) (Mg-mxa) , 
a-yjD^ = - {mx^+Z^) W , 

b^^D^ = (ly-M^) Z^+ (mXs-^Z^) , 

b^^D, = (Jy-W^)Zj^+(inXe+Z^)Wj^ , 

a2il>v = (m-Zj (mXs+Mji Z^ . 

a^zD^ = (in-Z^) {Mg-mx^) + (mx^+M^) (yn+Z,) , 

a23-^v = “ (m-Z^) W , 

b2^D^ = (m-Z^) Mf,+ (xoxq^M^) Zf,^ , 

^22^v ~ {HIXq'^’M^) Zj^ , 

d„ (w,g)D^ = ( ly-M^) I^+ {mxg+z^) I, , 
dg{w,g)Dy= (m-zj Ig+imXg+M^) I^, , 
c^{w,g)D^ = ily-M^) ntz^g^- (.mX(;+Z^) mz^wg , 

Cj (w,g)Dy = - (;n-Z^) mZgWg* (mx^+M^) mz^g^ . 


(2.9) 


In equations (2.5) through (2.8), the vehicle is assumed to be 
neutrally buoyant (W * B) , level (Xq * Xg) , and statically 
stable (zq > Zg) . The terms and Ig represent the cross flow 
drag integrals in equations (2.1) and (2.2), and ZQg Zq - Zg 
is the metacentric height. We can assume Zg to be zero, 
therefore z^g = Zq. 

The depth control system uses the linearized form of 
equations (2.5) through (2.8) with the linearization at the 
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nominal conditions of level attitude, ordered depth and 
forward speed. This gives the state space form of the depth 
control equations as 


6 = g , 


( 2 . 10 ) 


w ^ , ( 2 . 11 ) 


q ^ a^^U^w^a^^U^q^a^^ZsB^^b^U^^b . (2.12) 


i = -U^e*w , (2.13) 

where Uq is the nominal speed for gain selection, and the 
control inputs are defined as 

6 , = 6 , 

(2.14) 

i?2 ~ •^21'*'*^22 ' 

where a is the control surface coordination gain. This 
produces a linear full state feedback control law of 

6 = k^Q+k2W+k^q+k^z , (2.15) 

where the gains k^, kj, k 3 , and k^ are computed to give the 
closed loop system, equations (2.10) through (2.13), the 
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desired dynamics. Since the desired characteristic equation 
has the general form, 


= 0 , (2.16) 

the controller gains are computed by equating coefficients of 
the actual and desired characteristic equations, 

(^ 2312 - 2 ?! 322 ) ^0^2* <i>1^21-^2^l0 
- “®2“^23^GB‘*‘^^11^22"^21^12^ ' 

(2.17) 

(■^2^11 ”•^1^21^ ^•^1^23”■^2^13^ ^ GB ^ 0^^2 

+ (2)2+jbi322--b23l2) ^^12^21-^23^11^ ^GB^O ' 

[ (i3j^32i~.b23j^j_) Uq** (•^ i ^ 23 ~-^ 2 ^ 213 ^ ^GB^O^^ ~ ®0 ' 

[ref. 7] 
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Figure (2.1) - Geometrical representations of the basic 
definitions for the equations of motion. 
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III. STABILITY 


A. BIFURCATION ANALYSIS 

In system dynamics, the classical definition of stability 
states, that the real part of all the eigenvalues of the 
system must be negative. Therefore, our initial 
investigations into the stability of the SUBOFF Model was to 
find those eigenvalues whose real parts cross the imaginary 
axis. We used the bifurcation analysis program, included as 
Appendix A, to calculate the eigenvalues of the system. 

By linearizing the equations of motion, equations (2.1 
through 2.4), the state space equations of the dynamic system 
can be written in the form 


where 


u=-JCr 


(3.2) 


and K is the matrix of controller gains, as calculated by pole 
placement in equation (2.17). The eigenvalues of the system 
are found by solving 

(3 3) 

detlA-S/(r-sI)=0. ' 
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In the bifurcation analysis program a pseudo-root locus 
method is employed where the time constant, T^,, is fixed. The 
constant T^. fixes to placement of the system poles at a given 
nominal forward speed Ug and then the model speed, U, is 
varied incrementally with the system eigenvalues calculated at 
each speed increment. When the real part of an eigenvalue 
changes sign between the limits of a speed increment a 
bisection method is employed to find the speed where the real 
part of the eigenvalue equals 0. 

For each point where the real part of an eigenvalue 
crosses the imaginary axis the associated T^, and U are plotted 
on a bifurcation map. This map delineates the regions of 
classical stability (all eigenvalues on the left hand half- 
plane) from the regions of instability (at least one 
eigenvalue in the right hand half-plane). A family of 
bifurcation maps were generated by varying nominal speed, Ug, 
initial sted^ility, Zgg, and control surface gain, a. 

B. TYPICAL RESXTLTS 

Figure (3.1) shows a typical bifurcation map with its five 
distinct regions. Region I is the area of classical 
stability. In region II there is one real positive eigenvalue 
which is indicative of a pitchfork bifurcation. Pitchfork 
bifurcations of this model were previously examined by Reidel 
[ref. 7]. Regions III, IV, and V have at least one pair of 
complex conjugate eigenvalues with a positive real part. This 
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would indicate that there should be an unstable oscillatory 
behavior for the model. 

C. SIMULATIONS 

An extensive set of simulations were run to verify the 
bifurcation map's prediction of system stability. While the 
results of simulations showed the demarcation between the 
stable and unstable regions, the simulations found that the 
linear bifurcation analysis failed to predict the method of 
departure from stability- Five points (a through e) as shown 
on fig. (3.1) were chosen to illustrate the model's behavior 
in the regions of interest- Tc±>le (3.1) lists the eigenvalues 
found at each of these points and Table (3.2) shows the 
eigenvalues associated with the exact bifurcation point near 
points (b through d). Note that the eigenvalues are given in 
dimensional terms while all other information in the tables 
are non-dimensionalized. 

Point a is in the region of stability and fig (3.2) shows 
a rapid convergence to nominal stability. Simulations run at 
points b, c, and d are shown on figures (3.3), (3.4), and 
(3.5). These points are spaced less than 2.2 kts apart but 
they have a huge difference in their bounded oscillatory 
behavior. Table (3.1) also lists the measured and theoretical 
periods for these three points. The theoretical period is 
computed by dividing 2w by the imaginary part of the 
eigenvalue with the largest real part, and then non- 
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dimensionalize it by the ratio U/L, where U = speed of the 
vehicle (in feet per second) at each point, and L = vehicle 
length (13.9792 ft). An interesting point to note is the 
behavior at point d, where the measured period is 
approximately one half of the theoretical period. The 
dominant period at 80 nd sec is associated with the creation 
of the limit cycle at the bifurcation point while the 155.4 nd 
sec period is part of the chaotic response found at point d. 
Table (3.2) shows the eigenvalues, measured, and theoretical 
periods at the exact bifurcation points. Finally a simulation 
run at point is shown in figure (3.6) and demonstrates an 
unbounded departure from stability. 

It is evident that a more detailed analysis had to be 
performed, and that a Hopf bifurcation analysis should be 
used. 






TIME CONSTANT (non-dimensional) 


BIFURCATION MAP FOR UO = 9.0 FPS 



0 0.5 1 1.5 2 2.5 3 

VELOCITY (non-dimensional) 


Figure (3.1) - A typical bifurcation map showing the five 
distinct regions. 
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Theta (decrees) 


COMPARATIVE SIMUUTION PLOT 



Time (noodimensional) 


Figure (3.2) - An example of the steLble response found in 
region I. This corresponds to Point a in fig. (3.1). 
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Theta (degrees) 


COMPARATIVE SIMUUTION PLOT 



Figure (3.3) - This figure corresponds to Point b in fig 
(3.1) euid shows the oscillatory behavior in region III. 
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Tkela (degrees) 


COMPARATIVE SIMUUTION PLOT 



Figure (3.4) - This simulation corresponds to Point c in 
fig. (3.1). Note the distinct difference in the nature of 
oscillations between this point and Point b. 
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COMPARATIVE SIMULATION PLOT 



Figure (3.5) - This simulation corresponds to Point d in 
fig. (3.1). Note the magnitude of the oscillations associated 
with region IV. 
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COMPARATIVE SIMUUTION PLOT 



Time (noodimensionel) 


Figure (3.6) - The simulation run at Point e in region V 
of fig. (3.1) shows the unbounded behavior found in this 
region. 
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Table (3.1) - A listing of the points in fig. (3.1) and 
their associated eigenvalues, and periods. 
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1 Points 

Eigenvalues 

U 

(fps) 

Measured T 
(nd sec) 

Theoretical T 
(nd sec) 

b 

-0.6070, -0.0002 
0.0000 ± 0.4506 

O. 6 IU 0 

6.11 

5.57 

c 

-0.5710, -0.0054 
0.0000 ± 0.2884 

0.93Uo 

13.13 

13.18 

d 

-0.2675 ± 0.2063 
0.0000 ± 0.0558 

I.O 2 U 0 

74.93 

74.67 


Table (3.2) - Listing of the eigenvalues and periods of 
the exact bifurcation points associated with the points in 
Fig. (3.1). 
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IV. HOFF BIFURCATION 


A. INTRODUCTION 

By definition a Hopf bifurcation occurs when a pair of 
con^lex conjugate eigenvalues cross into the right hand half- 
plane. When this occurs the system will deviate from a steady 
state solution in an oscillatory manner. This deviation can 
be either supercritical or subcritical. 

For the supercritical case shown in fig. (4.1), stable 
limit cycles form after straight line stability is lost. 
Assume that a parameter, w(t), is varying. When t is less 
than t„it of eigenvalues of the system are located in 
the left hand half-plane and the system is nominally stable. 
At t equal to t^rit' ® complex conjugate pair moves into the 
right hand half-plane and forms the stable limit c/cle. As 
the distance w(t) increases from w(t 5 ,rit) / where the difference 
D = w(t) - w(t(,rit) ' amplitude of the limit cycle will also 
increase. If D remains small then the system will remain near 
the nominal steady state solution. 

Fig. (4.2) shows the subcritical case with unstable limit 
cycles being generated prior to the critical point being 
reached. Thus, as w(t) approaches w{t„it) f^® system could 
deviate from the nominal steady state solution and converge to 
a large amplitude limit cycle before the nominal system loses 
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stcUoility. Hence a random disturbance can cause a nominally 
stable system, one with all of the eigenvalues in the left 
hand half-plane, to exhibit oscillatory behavior. Once w(t) 
equals the nominal system becomes unstedDle and a 
discontinuous increase in the amplitude of oscillation is 
seen. This is because there are no nearby stable attractors 
for the system to converge to. 

A system design must make a distinction between these two 
types of bifurcations because of the disparate nature of the 
losses in stability. Thus the designer cannot rely on a 
linear approximation and must use higher order approximations 
(3^'^ order) of the equations of motion to adequately analyze 
his/her system. 

B. THIRD ORDER APPROXIMATIONS 

The nonlinear equations of motion are expanded using a 
third order Taylor series approximation near the nominal 
steady state, x = [0]. The control law is then modelled as, 

e' . S„.tanh(.^) , 


where is the saturation angle of the control plane input. 
Using the same approximation for the control law as the 
equation of motion, 6' becomes 


22 




1 


(4.2) 


6' = fi- 




Therefore, the state ecfuations can now be written as 

± = A* + gr(x) , (4.3) 


where 


x= [6 w g 2 ]“^ , 


(4.4) 


and the higher order terms are, 

^ 1 = 0 , 


(4.5a) 


^2 = bit;263(0,w,g,z) - -laijZg^e^ 

® (4.5b) 

+ c^^g^ + Cjjwg , 


^3 = Jb2U^dj(e,ur,g,z} - 4 -^z3^gb^^ 

O 

+ Cjiwg + c^zg^ , 


(4.5c) 


s. - 4'^' * • 


(4.5d) 
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The term 63 contains the third order terms derives from 
sxibstituting 6 into 5' . 

Defining T as the matrix of eigenvectors of A evaluated at 
the Hopf bifurcation point, then the transformation, 

X = Tz , (4.6) 


transforms the system into a canonical form 

i = T-^ATZ + T-^g(Tz) . (4.7) 


At the bifurcation point 


T'^AT 


0 -<i)o 0 0 

(Oq 0 0 0 

0 0 p 0 

0 0 0 g 


(4.8) 


with wq > 0 and p,q < 0. Z3 and correspond to p and q and 
are asymptotically stable. Center manifold theory states that 
the stcd 5 le coordinates Z 3 , z^ can be expressed as polynomials 
in the critical coordinates z^, Z 2 and this relationship is at 
least second order. Therefore, we can write, 

Zj = * “ 3 ^ 2 ^ ' (4.9a) 


and 


^4 = PlV P2^1^2 P3^2^ • 


(4.9b) 
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The coefficients and can be computed as follows: We 
differentiate ec[uations (4.9) with respect to z, 

+ P2(^1^2 ^ ^P3^2'^2 ' 

and substitute ±-^ » '‘^ 0^2 • *2 “ “ 0^1 • Therefore, 

(Sffj - 2oti) WqZ^Zj - ttjWoZj^ , (4.10a) 


and 


= Pa^oV + <2P3 - 2pi)&)oZiZ2 - P2 «o^2 ^ ■ {4.10b) 


The third and fourth equations of (4.7) are written as, 


3 _ |P 0| ^3 
■4 " lo <31 ^4 


+ ID] , 


(4.11) 


where, 

[D] = T-^g^iTz) , 

and 92 (Tz) contains the second order terms of g(Tz). We 
substitute equations (4.9) into (4.11) and equate coefficients 
with (4.10) . In this way we get a linear system in a^, and 
jS^. From this we can write the two dimensional state space 
equations as 
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{4.12a) 


Zi = -0)0^2 + + -^^13^1^2" ^ -^14^2' 

+ Pll^l^ + Pi 2^1^2 Pl3^2^ ' 


and 

Zj = (OqZi + r2iZ,3 + + r23ZiZ22 + 

(4.12b) 

■*■ Pai^l^ '*' P22'^l'^2 '*’ P23'^2^ ' 

where the coefficients r^j and p^j are derived from equations 
(4.8) . 

These equations are only valid exactly at the Hopf 
bifurcation point. For speeds U in a region near the 
bifurcation point these equations become 

Z 3 = a'ezj - (Wq + oa'e) Z 2 + F^iz^.z^) , (4.13a) 


and 


±2 = (<*>o ■*■ «'€)z 3 + a'czj + F 2 (z^,Z 2 ) , (4.13b) 


where: a', «' are the derivatives with respect to U of the 
real and imaginary parts of the critical complex conjugate 
pair of eigenvalues evaluated at the bifurcation point: e is 
the difference in U from the critical value; and, the 
nonlinear functions and F 2 are 
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(4.14a) 



Zj = i2sin0 (4.15b) 

equations (4.13) become 

R - a'eR + Fi(i?,0)cos© + F2(J?,0)sin0 (4.16a) 

and 

i30 = (<i) + «'€)/? + F2(i?,0)cos0 - Fi(i?,0) sin© . (4.16b) 

Equation (4.16a) then yields 

R = o'cF + P(©)F^ 0(0) . (4.17) 

By averaging equation (4.17) over one cycle we can obtain an 
equation with constant coefficients. Defining 
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(4.18) 


K = ^r^*P(e)d9 , 

2n Jo 

and 

L = -^f^*C)(e)c» , (4.19) 

2-nJo 

and carrying out the integration we obtain 

L = 0 , (4.20) 

and 

K = — + rj3 + ^22 • (4.21) 

which reduces equation (4.17) to 

R = a'cJ? + KR^ . (4.22) 

The existence and stability of the limit cycle is 
determined by analyzing the equilibrium points of the averaged 
equation (4.22), which correspond to periodic solutions in z-^ 
Z 2 as seen in the coordinate transformation equations (4.15) . 
From equation (4.22) we can see that two conditions exist, 

1) If a' > 0, then: 

(a) If K > 0 then unstable limit cycles coexist with the 
stable equilibrium for € < 0. 
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(b) If K < 0 then stable limit cycles coexist with the 
iinsteUDle equilibrium for e > 0. 

2) If ex' < 0, then: 

(a) If K > 0 the unstable limit cycles coexist with the 
stable ecjuilibrixim for e > 0. 

(b) If K < 0 then stcUale limit cycles coexist with the 
unstcUDle equilibrium for e <0. 

From this criteria, by confuting the nonlinear coefficient 

K we can use it to distinguish the two different types of Hopf 

bifurcations: 

• Supercritical if K < 0; 

• Subcritical if K > 0. 

[ref. 6] 

C. RESULTS 

A typical bifurcation map of stability for the SUBOFF 
Model is shown in fig. (4.3). This map is characterized by 
the Pitchfork bifurcation curve (P) and the three Hopf 
bifurcation curves (HI, H2, and H3) . The nature of curve P 
was previously analyzed and those results (Riedel, 1993) are 
reconfirmed in this study. 

Curve HI is characterized by a weak supercritical branch 
(a -* b) at low nominal speeds, Uq. As Uq increases this 
branch develops a weak to moderate subcritical behavior with 
K between 0 and 10^. The second branch of Hi (b -*■ c) has a 
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consistent moderate subcritical behavior with K on the order 
of 10 ^. 

Cusp (C) marks the intersection of curve H2 with curve H3. 
The cusp is highly dependent on both Uq and initial stability 
Zqb- 

As Zgg, for a given Uq, increases curve H2 (d -♦ e) shifts 
from a very weak subcritical nature with K between +10'^ and 
1 to a very weak supercritical nature with K between -l and - 
10'^. With a lower Ug and/or higher Zq^, point e moves down 
in the time constant and may not intersect curve H3. 

Curve H3 (f -* g) is a strong supercritical bifurcation 
with K values between -10^ and -10^. The position of H3 is 
independent of Ug# initial stability, and control surface 
coordination gain. 

Finally, curve H4 is a strong subcritical branch with K 
having values between 10^ and 10®, Because of this highly 
subcritical behavior, H4 can dominate and obscure the stable 
region at speeds greater than U/Ug = 1. 

Figure (4.4) plots the K values for a representative 
bifurcation map shown in fig (3.1) . Note the predicted super- 
and subcritical branches associated with fig (4.3) . Point (a) 
is inside the stcdjle region I and the numerical simulations 
seen in fig. (3.2) converge to zero. Point (b) is located in 
the unstable region, immediately after a supercritical 
bifurcation. As a result, small amplitude limit cycle 
oscillations have developed. The same is true as we move 
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towards point (c) although we expect a faunily of unstable 
limit cycles around this point as a result of the subcritical 
bifurcation. This point is further explored in the next 
section. As we approach point (d) a family of stable limit 
cycles is generated but its behavior is influenced by the 
previously developed unstadDle limit cycles. The real part of 
the critical pair of eigenvalues is becoming positive and 
relatively large, see fig. (3.1), which means that the 
cunplitudes of these stcdJle limit cycles are expected to be 
significantly higher, a result which is confirmed by fig. 
(3.5). The imaginary parts of the critical pair of 
eigenvalues are also changing very fast in this region. This 
means that the period of these limit cycles must be computed 
based on the value of the imaginary part right at the 
bifurcation point, rather than its value at the specific 
parameter point. This was observed previously in tables (3.1) 
and (3.2). Point (e) is in the strongly subcritical region V, 
thus we see the rapid divergence from stability as seen in 
fig. (3.6). 

D. SIMULATIONS 

The response of the submersible was simulated using an 
Adams-Bashforth integration scheme in Fortran program coding 
and is included in Appendix C. The contro? law (2.15) and 
gains (2.17) discussed in Chapters II and III were used. Non- 
dimensional ships speed and control system time constant. 
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nominal speed, initial stability, and control surface 
coordination gain were the input values to the system. A 
nominal 0.1 ft/sec vertical speed was used as the external 
initial disturbance. 

The simulations were used to compare the Hopf bifurcation 
data in two ways: 

1) By confirming sub-/ supercritical behavior predicted by 
the K factor, and; 

2) By comparing the predicted to simulated period of 
oscillations. 

In fig. (4.5) we have plotted the stable equilibria and 
the stable and unstable limit cycles from our example first 
used in fig. (3.1). Figure (4.5) clearly shows the predicted 
sub- and supercritical behavior that the K values predicted. 
The important features to note in fig. (4.5) are: 

1. Unstable limit cycles found with the subcritical Hopf 
bifurcation; 

2. The divergence of tne amplitudes as the velocity moves 
away from the bifurcation point velocity, and; 

3. The rapid divergence of the right most bifurcation and 
its quick and abrupt termination. 

To see why the abrupt termination occurs we must look at 
the root locus plot as shown in fig. (4.6) . The parameter in 
this the third bifurcation terminates when there is a break in 
point and the complex conjugate pair, which had given the 
oscillatory response, becomes two positive real eigenvalues. 
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These two eigenvalues now make the system completely unstable. 
The splitting of this complex conjugate pair into two real and 
positive eigenvalues can be better seen in fig (4.7). 
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SUPERCRITICAL HOPF BIFURCATION 



Figure (4.1) - An example of a supercritical Hopf 

bifurcation where the system has a nominal steady state 
solution until point C, which occurs at After t^^,. the 

system becomes unstable but forms a staJale limit cycle. 


34 



SUBCRITICAL HOPF BIFURCATIONf 



Figure (4,2) - This figure shows a subcritical Hopf 

bifurcation where the system loses stcUaility prior to reaching 
point C, which again occurs at t^rif 
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A TYPICAL BIFURCATION MAP 



Figure (4.3) - An annotace<3 bifurcation map for the SUBOFF 
Model. 


1 
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BIFURCATION MAP FOR UO = 9.0 FPS 




a 'b 
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+ + 


Zg = 0.4, Alpha = 0.0 


IV Pitchfork bifurcation (-) 


K value < 0 (.) 


K value 0 { + + 4. + + ) 
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Figure (4.4) - A plot of the K values associated with fig. 
(3.1). 
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UMIT CYCLES FOR UO = 9.0.Zg = 0.4. Alpha = 0.0 



Velocity (non-dimensional) 


Figure (4.5) - This shows the amplitude response for a 
speed range encompassing the Hopf bifurcation points that are 
shown in fig. (3.1). 
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associated with fig. (4.5). 




Real Axis 


REAL PART OF THE EIGENVALUES VS. U 
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0.2 
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Velocitj (non-dimensional) 

Figure (4.7) - A plot of the real parts of the eigenvalues 
plotted in fig. (4.6). 
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V. APPLICATIONS 


A. CONTROL SYSTEM PARAMETERS 

From the typical bifurcation maps we can see that a region 
of stability is created between the pitchfork and Hopf 
bifurcation boundaries. For the control system designer the 
limits or parameters must be defined prior to starting the 
design. By maximizing the region of stability we can give the 
designer the most leeway in his/her work. There are three 
parcuneters that we can use to change the bifurcation maps, 
nominal speed, initial stability, and the control surface gain 
coefficient. 

First we will look at changing the nominal speed. In fig. 
(5.1) we have plotted three curves for nominal speeds of 3.0, 
9.0, and 15.0 fps. We can see that although the pitchfork 
line moves to the left in dimensional speeds this line remains 
nearly constant with a dimensional stern planes reversal 
occurring at 1.2 kts. The high speed Hopf boundaries (U/Uq > 
1.0) move apart as the nominal speed increases. The 
effectiveness of increasing Uq is limited in the upper branch 
by the fixed position of the H3 line with the maximum 
practical T^, achieved at a Uq = 9.0 fps. In the lower arm 
there is no increase in the stability area after Uq * 9.0 fps 
therefore any increase in Uq is pointless. For the low speed 
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Hopf curves (U/Ug < 1.0) we quickly lose our margin of 
stability as Ug increases thus necessitating further changes 
to regain the lost area of stability. 

The next parameter change we examined was varying the 
initial stability. Figure (5.2) shows the effect of 
increasing Zg from 0.1 to 0.4 ft. The subcritical H4 branch 
remains constant while the upper high speed Hopf branch moves 
down effectively decreasing the area of stability. The low 
speed Hopf curves move up to increase the low speed area of 
stability. We can see that the additional loss in area is by 
the movement of the pitchfork line to the right. At a Zg * 
0.4 ft stern planes reversal occurs at a dimensional speed of 
2.4 kts which is well within the currently accepted range of 
1.0 to 3.0 kts for modern submarines. Therefore in this 
model we would want to balance the initial staUaility to 
maximize the low and high speed areas. 

An increase in the control surface gain coefficient, a, is 
shown in fig. (5.3). Note that the low and high speed Hopf 
curves all move up in time constant. While the low speed Hopf 
curves give a large increase in stability, the high speed 
curves move up proportionally and there is no increase in 
stability area. This does allow the designer to shift the 
range of stable time constants without a loss of high speed 
stability. The pitchfork line will move to the left until it 
equals zero when a = 1.0. 
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In order to examine what happens at the extremes of the 
design options we can look at the low nominal speed (3 fps) 
bifurcation maps. Figure (5.4) shows the typical bifurcation 
map such as the one we have previously discussed. As the 
metacentric height is increased there are significant changes 
in the nature of the bifurcation curves- In fig. (5.5) we see 
that the pitchfork bifurcation line has moved significantly to 
the right and has intersected the low speed Hopf bifurcation 
curve. This intersection along with the merger of the H2 and 
H4 curves have combined to reduce the region of stability to 
a negligible portion of the map. 

A further increase in the metacentric height as shown in 
fig. (5.6) demonstrates a dreunatic change in the nature of the 
stability of the model. The low speed region has two 
hyperbolic-like Hopf bifurcation curves (the upper curve 
occurs well above the region of interest) bounding the lower 
and upper limits of stability. For the speeds U/Uq > 1.0, the 
pitchfork bifurcation line now intersects the H2 curve and 
has changed from a supercritical to a subcritical pitchfork. 

Figure (5.7) is an example of the effect just described 
and how it can occur at higher nominal speeds. This shows 
that although initial stability is necessary for overall 
stability, if the metacentric height becomes too large it can 
have an adverse affect on the performance of the submarine. 
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B. SUBNBRIIIB DESIGN EVALUATION 

For all of the simulations up to this point the moment and 
inertia coefficients of the control surfaces have been .5x of 
those listed in the SUBOFF Model report. By using the reduce 
effect control surface input we have been able to show stable 
simulations in all of the mapped stable regions of the 
bifurcation maps. 

Linear bifurcation methods fail to predict the change in 
system response for changes in control surface coefficients. 
The bifurcation maps for .lx to l.Ox the control surface 


coefficients 

are exactly 

the same, 

in other 

words 

the 

bifurcation 

points are 

independent 

of the 

size 

and 


effectiveness of the control surfaces. Therefore we must 
examine the K values in order to predict the response of the 
model. 

Figure (5.8) shows the change in stability for the model 
with and increase in the control surface coefficients. The 
area of lost stability is indicated by the shaded portion of 
the map. This loss of stability is caused by the shift of the 
H2 curve from weak to moderately supercritical to a strongly 
subcritical curve. With two strongly subcritical curves in 
the high speed region (U/Uq > 1.0) the possibility of 
subcritical capture is greatly increased. This effect is 
confirmed by running an extensive set of simulations and 
mapping the change in stable to unstable response. We must 
note that this instability occurs in a region that has four 
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eigenvalues with negative real parts where linear control 
system design would not predict an instaUoility. 

With this method a design can be tested and then modified 
using the design goals as determined by the nonlinear response 
in the simulations. This can change a completely unstable 
model such as the SUBOFF Model into one where a large margin 
of stcdiility and control system latitude is design into it. 
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Time Constaat (noa-dimenslonal) 


COMPARATIVE BIFURCATION MAP: INCREASING Zg 



Figure (5.2) - The effects of changing initial stability, 
Zg, on the bifurcation maps. 
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Time Constant (non-dimensional) 


COMPARATIVE BIFURCATION MAP: INCREASING Alpha 



Figure (5.3) - The effects of changing the control surface 
gain coefficient, a, on the bifurcation maps. 





TIME CONSTANT (non-dimensional) 


STABIUTY/BIFURCATION MAP FOR UO = 3 FPS 



Figure (5.4) - A low nominal speed bifurcation map. 


49 



TIME 



0 0.5 1 1.5 2 

VELOCITY (non-dimensional) 


Figure (5.5) - This demonstrates the effects of an 

increase in the metacentric height for low nominal speeds. 
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BIFURCATION MAP FOR UO = 3.0 FPS 


^^Stable Region 

Ox 


VELOCITY (non-dimensional) 


Figure (5.6) - As the metacentric height increases it has 
adverse effect on the stability of the siobmarine. 










TIME CONSTANT (non-dimensional) 


bifurcation map for UO s 8.0 FPS 



Figure (5.7) - An example of increasing the metacentric 

height for a higher nominal speed model. 
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Figure (5.8) - 
moment and inertia 
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VI. CONCLUSIONS AND RBCG9MKENDATIONS 


A. CONCLUSIONS 

The application of a Hopf bifurcation analysis to a 
submarine design can be an effective tool in the design 
evaluation and modification phase. These methods when paired 
with programs that generate hydrodynamic coefficients for a 
sxibmarine will save time, effort, and money by reducing the 
amount of model testing necessary to validate a design. An 
effective set of control system design parameters will be 
generated in the process that will be optimal for the final 
design of a submersible. 

Secondly, this system will give the limits of the range of 
metacentric heights that will maintain steQsility for the full 
range of speeds of the design. As was shown changes in the 
metacentric height can have a dramatic affect on stability. 

Finally, an evaluation of the need for bow planes or 
forward control surfaces are actually to maintain control of 
the submarine. If the forward planes can be eliminated a 
potential source of noise would also be eliminated along with 
simplifying the design of the control system. 
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B. RBCQIOISMDATIQITS 


First, modify the programs to evaluate the performance of 
the svibmarine for the effects of external forces such as, wave 
effects, currents, and free surface effects. 

Second, perform a systematic study of the bifurcations at 
conditions other than nominal in trim conditions, such as out 
of neutral buoyancy (heavy/light), or out-of-trim fore and aft 
conditions. 
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AFPBMDIZ A - BIFORCATION ANALYSIS PROGRAM 


PROGRAM BIFURl.FOR 

BIFURCATION ANALYSIS FOR THE DARPA SUBOFF MODEL 

PARAMETERS ARE: TC VS. U 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2,K3,K4,L,MQDOT, 

& MWDOT,MQ,MW,MDS,MDB,MD, 

& MASS,IY,WF,pi 

DIMENSION A(4,4).FV1{4),IV1(4),ZZZ(4,4),WR(4),WI(4) 

OPEN (11,FILE-'BIFl.RES',STATUS-'NEW') 

OPEN (12,FILE-'BIF2.RES',STATUS-'NEW') 

OPEN (13,FILE-'BIF3.RES',STATUS-'NEW') 

WEIGHT-1556.2363 


BUO 

-1556.2363 

L 

-13.9792 

lY 

-561.32 

G 

-32.2 

MASS 

-WEIGHT/G 

RHO 

-1.94 

XB 

-0.0 

ZB 

-0.0 

pi 

-3.14159 


WRITE 

(*,*) ' 

ENTER MIN, MAX, AND # OF INCREMENTS IN TIME 
CONSTANT' 

READ 

(*,*) 

TCMIN,TCMAX,ITC 

WRITE 

(*,*) 

'ENTER MIN, MAX, AND # OF INCREMENTS IN 
VEHICLE SPEED' 

READ 

(*,*) 

UMIN,UMAX,IU 

WRITE 

(*,*) 

'ENTER NOMINAL SPEED' 

READ 

(*, *) 

UO 

WRITE 

(*,*) 

'ENTER XG AND ZG' 

READ 

(*,*) 

XG,ZG 

TCMIN 

= 0.1 


TCMAX 

- 9.0 


ITC 

= 200 


lU 

= 200 


UO 

H 

I-* 

to 

• 

o 


XG 

- 0.0 


WRITE 

(*,*) 

'ENTER MIN, MAX, OF VEHICLE SPEED' 

READ 

(*,*) 

UMIN,UMAX 

WRITE 

(*,*) 

'ENTER ZG' 

READ 

(*,*) 

ZG 
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c 

ZGB-Z6-ZB 
TCMIN-TCMIN* L/UO 
TCMAX-TCMAX* L/UO 
UMIN «UMIN*UO 
UMAX -UMAX*UO 
C 

ZQDOT--6.3300E-04*0.5*RHO*L**4 
ZWDOT--1.4529E-02*0.5*RHO*L**3 
ZQ = 7.5450E-03*0.5*RHO*L**3 
ZW --1.3910E-02*0.5*RHO*L**2 
ZDS «0.5*(-5.6030E-03*0.5*RHO*L**2) 

ZDB -0.5*(-5.6030E-03*0.5*RHO*L**2) 

MQDOT--8.8000E-04*0.5*RH0*L**5 
MWD0T--5.6100E-04*0.5*RHO*L**4 
MQ --3.7020E-03*0.5*RHO*L**4 
MW - 1.0324E-02*0.5*RHO*L**3 
MDS -0.5*(-2.4090E-03*0.5*RHO*L**3) 

MDB -0.5* ( 2.4090E-03*0.5*RHO*L**3) 

C 

ALPHA--0.0 
ZD-ZDS +ALPHA* ZDB 
MD-MDS +ALPHA*MDB 
C 

DV =(MASS-ZWDOT)*(lY-MQDOT) 

& -(MASS*XG+ZQDOT)*(MASS*XG+MWDOT) 

AllDV-(lY-MQDOT)*ZW+(MASS*XG+ZQDOT)*MW 

A12DV- (lY-MQDOT) * (MASS+ZQ) + (MASS*XG+ZQDOT) * (MQ-MASS*XG) 

A13DV--(MASS*XG+ZQDOT)*WEIGHT 

BIDV =(lY-MQDOT)*ZD+(MASS*XG+ZQDOT)*MD 

A21DV- (MASS-ZWDOT) *MW+ (MASS*XG+MWDOT) *ZW 

A22DV-(MASS-ZWDOT)*(MQ-MASS*XG) 

& +(MASS*XG+MWDOT)*(MASS+ZQ) 

A23DV--(MASS-ZWDOT)*WEIGHT 
B2DV =(MASS-ZWDOT)*MD+(MASS*XG+MWDOT)*ZD 
C 

All-AllDV/DV 
A12-A12DV/DV 
A13-A13DV/DV 
A21-A21DV/DV 
A22-A22DV/DV 
A23-A23DV/DV 
B1 -BIDV /DV 
B2 -B2DV /DV 
C 

EPS -l.D-5 
ILMAX-1500 
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c 

DO 1 I-1,ITC 

WRITE (*,2001) I,ITC 

TC-TCMIN+(I-l)*(TCMAX-TCMIN)/(ITC-1) 

POLE-1.0/TC 
ALPHA3-4.0*POLE 
ALPHA2-6.0*POLE**2 
ALPHAl-4.0*POLE**3 
ALPHAO- POLE**4 
K4-ALPHA0/((B1*A21-B2*A11)*U0**4 
& +(B1*A23-B2*A13)*ZGB*U0**2) 

A2M-B1*U0**2 

A3M-B2*U0**2 

A0M-- (A11+A22) *U0-ALPHA3 
B1M»B2*U0**2 

B2M-(B2*A12-B1*A22)*U0**3 
B3M- (B1*A21-B2*A11)*U0**3 

BOM-(A11*A22-A21*A12)*U0**2-A23*ZGB-ALPHA2-B1*U0*U0*K4 
CIM-(B2*A11-B1*A21)*U0**3 
C2M- (B1*A23-B2*A13)*ZGB*U0**2 
COM- (A13*A21-A23*A11)*ZGB*U0+ALPHA1 
& -(B2+B1*A22-B2*A12)*K4*U0**3 

K2-C1M*B0M*A3M- B1M*C0M*A3M- C1M*B3M*A0M 
K2-K2/ (C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

Kl-(C0M-C2M*K2)/CIM 
K3- {A0M-A2M*K2) /A3M 

DO 2 J«1,IU 

U-UMIN+ (J-1) * (UMAX-UMIN) / (IU-1) 

A(1,1)=0.0D0 

A(1,2)-0.0D0 

A(1,3)-1.0D0 

A(1,4)-0.0D0 

A(2,1)-A13*ZGB+B1*U*U*K1 

A(2,2)»A11*U +B1*U*U*K2 

A(2,3)-A12*U +B1*U*U*K3 

A (2,4)= B1*U*U*K4 

A(3,1)-A23*ZGB+B2*U*U*K1 

A(3,2)-A21*U +B2*U*U*K2 

A(3,3)=A22*U +B2*U*U*K3 

A(3,4)= B2*U*U*K4 

A(4,l)=-U 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)-0.0D0 

<» 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 
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IF (J.GT.l) GO TO 10 

DEOSOO-OEOS 

UOO -U 

LL-0 

GO TO 2 

DEOSNN>DEOS 

UNN -U 

PR-DEOSNN*DEOSOO 

IF (PR.GT.O.DO) GO TO 3 

LL-LL+1 

IF (LL.GT.3) STOP 1000 

IL-0 

UO-UOO 

UN-UNN 

DEOSO-DEOSOO 

DEOSN-DEOSNN 

UL-UO 

UR-UN 

DEOSL>DEOSO 

DEOSR-DEOSN 

U»(UL+UR)/2.D0 

A(1,1)-0.0D0 

A(1,2)-0.0D0 

A(1,3)-1.0D0 

A(1,4)-0.0D0 

A(2,1)«A13*ZGB+B1*U*U*K1 
A(2,2)«A11*U +B1*U*U*K2 

A(2,3)-A12*U +B1*U*U*K3 

A(2,4)- B1*U*U*K4 

A(3,1)-A23*ZGB+B2*U*U*K1 
A(3,2)-A21*U +B2*U*U*K2 

A(3,3)=A22*U +B2*U*U*K3 

A(3,4)= B2*U*U*K4 

A(4,l)=-U 
A(4,2)=1.0D0 
A(4,3)=0.0D0 
A(4,4)«0.0D0 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 

DEOSM»DEOS 

UM-U 

PRL-DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 

UO»UL 

UN»UM 

DEOSO-DEOSL 

DEOSN=DEOSM 



IL-IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF-DABS(UL-UM) 

IF (DIF.GT.EPS) GO TO 6 

U-UM 

GO TO 4 

5 IF (PRR.GT.O.DO) STOP 3200 
UO-UM 
UN-UR 

DEOSO-DEOSM 

DEOSN-DEOSR 

IL-IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF-DABS(UM-UR) 

IF (DIF.GT.EPS) GO TO 6 
U-UM 

4 LLL-IO+LL 

WRITE (LLL,2002) U/U0,TC*U0/L 
3 UOO-UNN 

DEOSOO-DEOSNN 
2 CONTINUE 
1 CONTINUE 
C 

2001 ’'ORMAT (215) 

2002 FORMAT (4(3X,FIO.6)) 

END 

C 

SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4),WI(4) 

DEOS--1.0D+20 
DO 1 1-1,4 

IF (WR(I).LT.DEOS) GO TO 1 
DEOS-WR(I) 

IJ=I 

1 CONTINUE 
OMEGA-WI (IJ) 

OMEGA-DABS(OMEGA) 

RETURN 

END 
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APPENDIX B - HOPF BIFURCATION PROGRAM 


C 

C 

c 

c 


c 


c 


c 


PROGRAM HOPFIA.FOR 

EVALUATION OF HOPF BIFURCATION FORMULAE USING THE SUBOFF 
SUBMARINE MODEL 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L,IY,MASS,MQDOT,MWDOT, 

MQ,MW,MD,MDS,MDB,K1,K2,K3,K4, 

ALFAl,ALFA2,ALFA3, 

BETAl,BETA2,BETA3 

DOUBLE PRECISION M11,M12,M13,M14,M21,M22,M23,M24, 

M31 ,M32,M33,M34,M41,M42 , M43 ,M44, 

N11,N12,N13,N14,N21,N2 2,N2 3,N2 4, 
N31,N32,N33,N34,N41,N42,N43,N44, 

L21,L22,L23,L24,L31,L32,L33,L34, 

L41,L42,L43,L44,L25,L26,L27,L35, 

L36,L37,L21A,L22A,L23A,L24A,L31A, 
L32A,L33A,L34A 

DIMENSION A(4,4),T(4,4),TINV(4,4),FV1(4),IV1(4),YYY(4,4) 
DIMENSION WR(4),WI(4),TSAVE(4,4), 

1 TLUD(4,4),IVLUD(4),SVLUD(4) 

DIMENSION ASAVE(4,4),A1(4,4),A2(4,4) 

OPEN (10,FILE='HOPF.DAT',STATUS^'OLD') 

OPEN (20,FILE='HOPF.RES',STATUS^'NEW') 

WEIGHT=1556.2363 
lY = 561.32 
L = 13.9792 

RHO = 1.94 

G = 32.2 

XG = 0.0 

ZB = 0.0 

MASS =WEIGHT/G 
BOY ^WEIGHT 

ZQDOT =-6.3300E-04*0.5*RHO*L**4 
ZWDOT =-1.4529E-02*0.5*RHO*L**3 
ZQ = 7.5450E-03*0.5*RHO*L**3 
ZW =-1.3910E-02*0.5*RHO*L**2 
ZDS =0.5*(-5.6030E-03*0.5*RHO*L**2) 

ZDB =0.5*(-5.6030E-03*0.5*RHO*L**2) 

MQDOT «-8.8000E-04*0.5*RHO*L**5 
MWDOT =-5.6100E-04*0.5*RHO*L**4 
MQ «-3.7020E-03*0.5*RHO*L**4 
MW = 1.0324E-02*0.5*RHO*L**3 
MDS =0.5*(-2.4090E-03*0.5*RHO*L**3) 
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0.5*( 2.4090E-03*0.5*RHO*L**3) 


WRITE 

(*, 

1001) 


READ 

(*, 

*) 

ITOTAL 

WRITE 

(*, 

1002) 


READ 

(*, 

*) 

U0,ZG,DSAT 

WRITE 

(*. 

1004) 


READ 

(*, 

*) 

ZG 

WRITE 

(*. 

1003) 


READ 

(*, 

*) 

ALPHA 


Z6B»Z6-ZB 
UO =9.0 
DSAT = 0.4 

ZD=ZDS +ALPHA* ZDB 
MD-MDS+ALPHA*MDB 

DETERMINE [A] AND [B] COEFFICIENTS 


DV =(MASS-ZWDOT)*(lY-MQDOT) 

-(MASS*XG+ZQDOT)*(MASS*XG+MWDOT) 

A11DV=(lY-MQDOT)*ZW+{MASS*XG+ZQDOT)*MW 

A12DV= (lY-MQDOT) * (MASS+ZQ) + (MASS*XG+ZQDOT) * {MQ-MASS*XG) 

A13DV=-(MASS*XG+ZQDOT)*WEIGHT 

BIDV =(lY-MQDOT)*ZD+(MASS*XG+ZQDOT)*MD 

A21DV=(MASS-ZWDOT)*MW+(MASS*XG+MWDOT)*ZW 

A22DV=(MASS-ZWDOT)*(MQ-MASS*XG) 

+(MASS*XG+MWDOT)*(MASS+ZQ) 

A23DV=-(MASS-ZWDOT)*WEIGHT 

B2DV =(MASS-ZWDOT)*MD+(MASS*XG+MWDOT)*ZD 


A11=A11DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 
B1 =B1DV/DV 
B2 =B2DV/DV 

C11DV=(lY-MQDOT)*MASS*ZG 
C12DV=-(MASS*XG+ZQDOT)*MASS*ZG 
C21DV=-(MASS-ZWDOT)*MASS*ZG 
C22DV=(MASS*XG+MWDOT)*MASS*ZG 

C11=C11DV/DV 

C12=C12DV/DV 

C21=C21DV/DV 

C22=C22DV/DV 
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C 


DO 1 IT-1,IT0TAL 

WRITE (*,3001) IT,ITOTAL 
READ (10,*) U,TC 
TC-TC*L/U0 
U=U*U0 

CALCULATE THE GAINS 

POLE-l.O/TC 
ALPHA3-4.0*POLE 
ALPHA2-6.0*POLE**2 
ALPHAl-4.0*POLE**3 
ALPHAO- POLE**4 
K4-ALPHA0/((B1*A21-B2*A11)*U0**4 
& +(B1*A23-B2*A13)*ZGB*U0**2) 

A2M-B1*U0**2 

A3M=B2*U0**2 

A0M=-(A11+A22)*U0-ALPHAS 
B1M-B2*U0**2 

B2M- (B2 *A12 - B1*A22)*U0**3 
B3M-(B1*A21-B2*A11)*U0**3 

BOM- (A11*A22-A21*A12) *U0**2-A23*ZGB-ALPHA2-B1*U0*U0*K4 
CIM-(B2*A11-B1*A21)*U0**3 
C2M-(B1*A23-B2*A13)*ZGB*U0**2 
COM- (A13*A21-A23*A11) *ZGB*U0+ALPHA1 
& -(B2+B1*A22-B2*A12)*K4*U0**3 

K2=C1M*B0M*A3M-B1M*C0M*A3M-C1M*B3M*A0M 
K2=K2/ (C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

Kl=(C0M-C2M*K2)/CIM 
K3= (A0M-A2M*K2) /A3M 


EVALUATE NONLINEAR RUDDER EXPANSION COEFFICIENTS 

DTTW=-(l.D0/{3.D0*DSAT**2))*3.D0*K1*K1*K2 
DTWW=-(l.DO/(3.D0*DSAT**2))*3.D0*K1*K2*K2 
DTTQ=-(l.DO/(3.D0*DSAT**2))*3.D0*K1*K1*K3 
DTQQ--(l.DO/(3.D0*DSAT**2))*3.D0*K1*K3*K3 
DTTZ--(l.DO/(3.D0*DSAT**2))*3.D0*K1*K1*K4 
DTZZ--(l.DO/(3.D0*DSAT**2))*3.D0*K1*K4*K4 
DWWQ--(l.D0/(3.D0*DSAT**2))*3.D0*K2*K2*K3 
DWQQ--(l.DO/(3.D0*DSAT**2))*3.D0*K2*K3*K3 
DWWZ--(l.DO/(3.D0*DSAT**2))*3.D0*K4*K2*K2 
DWZZ--(l.DO/(3.D0*DSAT**2))*3.D0*K4*K4*K2 
DQQZ--(l.DO/(3.D0*DSAT**2))*3.D0*K4*K3*K3 
DQZZ--(l.D0/(3.D0*DSAT**2))*3.D0*K4*K4*K3 
DTWQ--(l.DO/(3.D0*DSAT**2))*6.D0*K1*K2*K3 
DTWZ--(l.DO/(3.D0*DSAT**2))*6.D0*K1*K4*K2 
DTQZ--(l.D0/(3.D0*DSAT**2))*6.D0*K1*K4*K3 
DWQZ--(l.D0/(3.D0*DSAT**2))*6.D0*K4*K2*K3 
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DTTT--(l.DO/(3.D0*DSAT**2))*1.D0*K1*K1*K1 
DWWW--(l.D0/(3.D0*DSAT**2))*1.D0*K2*K2*K2 
DQQQ--(l.DO/(3.D0*DSAT**2))*1.D0*K3*K3*K3 
DZZZ--(l.DO/(3.D0*DSAT**2))*1.DO*K4*K4*K4 

EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 

A(1,1)-0.0D0 

A(1,2)-O.ODO 

A(1,3)-1.0D0 

A(1,4)-0.0D0 

A(2,1)-A13*ZGB+B1*U*U*K1 
A(2,2)-A11*U +B1*U*U*K2 

A(2,3)-A12*U +B1*U*U*K3 

A(2,4)= B1*U*U*K4 

A(3,1)=A23*ZGB+B2*U*U*K1 
A(3,2)-A21*U +B2*U*U*K2 

A(3,3)=A22*U +B2*U*U*K3 

A(3,4)= B2*U*U*K4 

A(4,l)--U 
A(4,2)-1.0D0 
A(4,3)-0.0D0 
A(4,4)=0.0D0 
DO 11 1-1,4 
DO 12 J-1,4 

ASAVEd, J)-A(I, J) 

12 CONTINUE 
11 CONTINUE 

CALL RG(4,4,A,WR,WI,l,yyY,IVl,FVl,IERR) 

CALL DSOMEG(IEV,WR,WI,OMEGA,CHECK) 

OMEGAO-OMEGA 
DO 5 1=1,4 

T(i,i)= yyy(i,iEV) 

T(i,2)--yyy(i,iEV+i) 

5 CONTINUE 

IF (lEV.EQ.l) GO TO 13 

IF (IEV.EQ.2) GO TO 14 

IF (IEV.EQ.3) GO TO 15 

STOP 3004 

14 DO 6 1=1,4 

T(i,3)=yyy(i,i) 

T(i,4)=yyy(i,4) 

6 CONTINUE 
GO TO 17 

15 DO 7 1-1,4 

T(i,3)=yyy(i,i) 

T(i,4)=yyy(i,2) 

7 CONTINUE 
GO TO 17 

13 DO 16 1-1,4 

T(i,3)=yyy(i,3) 
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T(I,4)-YYY(I,4) 

16 CONTINUE 

17 CONTINUE 

CHECK THE EIGENVALUES 

WRITE(21,2003) TC*U0/L,WI(1),WI(2),WI(3),WI(4), 
& WR(1) ,WR(2) ,WR(3) ,WR(4) 

NORMALIZATION OF THE CRITICAL EIGENVECTOR 

INORM-1 

IF (INORM.NE.O) CALL NORMAL(T) 

INVERT TRANSFORMATION MATRIX 

DO 2 1-1,4 
DO 3 J-1,4 

TINVd, J)=0.0 
TSAVEd, J)=T(I, J) 

3 CONTINUE 
2 CONTINUE 

CALL DLUD(4,4,TSAVE,4,TLUD,IVLUD) 

DO 4 1-1,4 

IF (IVLUDd) .EQ.O) STOP 3003 

4 CONTINUE 

CALL DILU(4,4,TLUD,IVLUD,SVLUD) 

DO 8 1*1,4 
DO 9 J=l,4 

TINV(I,J)=TLUD(I,J) 

9 CONTINUE 

8 CONTINUE 

CHECK Inv(T)*A*T 

IMULT-1 

IF (IMULT.EQ.l) CALL MULT(TINV,ASAVE,T,A2) 

IF (IMULT.EQ.O) STOP 
P=A2 (3,3) 

Q=A2 (4,4) 

WRITE(21,*) P,Q 

DEFINITION OF Nij 

N11-TINV(1,1) 

N21-TINV(2,1) 

N31-TINV(3,1) 

N41-TINV(4,1) 

N12=TINV(1,2) 

N22=TINV(2,2) 

N32=TINV(3,2) 
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N42«TINV(4,2) 

N13-TINV(1,3) 

N23-TINV{2,3) 

N33-TINV{3,3) 

N43-TINV(4,3) 

N14-TINV(1,4) 

N24-TINV(2,4) 

N34-TINV{3,4) 

N44-TINV{4,4) 

C 

C DEFINITION OF Mij 

C 

M21-T(2,l) 

M31-T(3,l) 

M41=T(4,1) 

M12=T(1,2) 

M22=T(2,2) 

M32-T(3,2) 

M42-T(4,2) 

M13=T(1,3) 

M23=T{2,3) 

M33=T(3,3) 

M43=T(4,3) 

M14»T(1,4) 

M24=T(2,4) 

M34=T(3,4) 

M44=T(4,4) 

C 

C DEFINITION OF Lij 

C 

L25=C11*M31*M31+C12*M21*M31 

C 

L26=2*C11*M31*M32+C12*(M21*M32+M22*M31) 

C 

L27=C11*M32*M32+C12*M22*M32 

C 

L35=C22*M31*M31+321*M21*M31 

C 

L36=2*C22*M31*M32+C21*(M21*M32+M22*M31) 

C 

L37*=C22*M32*M32+C21*M22*M32 

C 

L21= DTTW*M11*M11*M21 + DTWW*M11*M21*M21 + 

& DTTQ*M11*M11*M31 + DTQQ*M11*M31*M31 + 

& DTTZ*M11*M11*M41 + DTZZ*M11*M41*M41 + 

& DWWQ*M21*M21*M31 + DWQQ*M21*M31*M31 + 

& DWWZ*M21*M21*M41 + DWZZ*M21*M41*M41 + 

& DQQZ*M31*M31*M41 + DQZZ*M31*M41*M41 + 

& DTWQ*M11*M21*M31 + DTWZ*M11*M21*M41 + 

& DTQZ*M11*M41*M31 + DWQZ*M21*M41*M31 + 
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& DTTT*M11*M11*M11 + DWWW*M21 *M21 *M21 + 

& DQQQ*M31*M31*M31 + DZZZ*M41*M41*M41 

TTW-Mll*Mll*M22+2.0*M11*M12*M21 
TWW«M12*M21*M21+2.0*M11*M21*M22 
TTQ-Mll*Mll*M32+2.0*M11*M12*M31 
TQQ-M12*M31*M31+2.0*M11*M31*M32 
TTZ-Mll*Mll*M42+2.0*M11*M12*M41 
TZZ-M41*M41*M12+2.0*M11*M41*M42 
WWQ-M21*M21*M32+2.0*M31*M21*M22 
WQQ-M22*M31*M31+2.0*M31*M32*M21 
WWZ-M21*M21*M42+2.0*M41*M21*M22 
WZZ«M2 2 *M41 *M41+2.0 *M41 *M4 2 *M21 
QQZ-M31*M31*M42+2.0*M41*M31*M32 
QZZ-M3 2 *M41*M41+2.0*M41 *M4 2 *M31 
TWQ-Ml 1 *M21 *M32 +M31 * (Ml 1 *M22 +M12 *M21) 
TWZ-M11*M21*M42+M41* (M11*M22+M12*M21) 
TQZ-Ml 1 *M41 *M3 2 +M31 * (Ml 1 *M4 2 +M12 *M41) 
WQZ=M21 *M41 *M32 +M31 * {M21 *M42 +M22 *M41) 

TTT»3.0*M11*M11*M12 
WWW«3.0*M21*M21*M22 
QQQ- 3.0 *M31 *M31 *M3 2 
ZZZ-3.0*M41*M41*M42 

L22=DTTW*TTW+DTWW*TWW+DTTQ*TTQ+DTQQ*TQQ+ 

& DTTZ*TTZ+DTZZ*TZZ+DWWQ*WWQ+DWQQ*WQQ+ 

& DWWZ*WWZ+DWZZ*WZZ+DQQZ*QQZ+DQZZ*QZZ+ 

& DTWQ*TWQ+DTWZ*TWZ+DTQZ*TQZ+DWQZ*WQZ+ 

& DTTT*TTT+DWWW*WWW+DQQQ*QQQ+DZZZ*ZZZ 

TTW«M12*Ml2*M21+2.0*M11*M12*M22 
TWW-Mll*M22*M22+2.0*M12*M21*M22 
TTQ-M12*M12*M31+2.0*M11*M12*M32 
TQQ-Mll*M32*M32+2.0*M12*M31*M32 
TTZ-M12*M12*M41+2.0*M11*M12*M42 
TZZ-Mll*M42*M42+2.0*M12*M41*M42 
WWQ-M22 *M22 *M31+2.0 *M21 *M22 *M32 
WQQ«M21 *M3 2 *M3 2+2.0 *M2 2 *M31 *M3 2 
WWZ-M22*M22*M41+2.0*M21*M22*M42 
WZZ»M21*M42*M42+2.0*M22*M41*M42 
QQZ-M32*M32*M41+2.0*M31*M32*M42 
QZZ-M31 *M4 2 *M4 2+2.0 *M3 2 *M41 *M4 2 
TWQ«M12*M22*M31+M32* (M11*M22+M12*M21) 
TWZ-M12*M22*M41+M42* (M11*M22+M12*M21) 
TQZ-Ml 2 *M4 2 *M31+M3 2 * (Ml 1 *M4 2 +M12 *M41) 
WQZ-M2 2*M42*M31+M32*(M21*M42 +M2 2 *M41) 

TTT-3.0*M11*M12*M12 
WWW-3.0*M21*M22*M22 
QQQ=3.0*M31*M32*M32 
ZZZ-3.0*M41*M42*M42 
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L23-DTTW*TTW+DTWW*TWW+DTTQ*TTQ+DTQQ*TQQ+ 

& DTTZ*TTZ+DTZZ*TZZ+DWWQ*WWQ+DWQQ*WQQ+ 

& DWWZ*WWZ+DWZZ*WZZ+DQQZ*QQZ+DQZZ*QZZ+ 

& DTWQ*TWQ+DTWZ*TWZ+DTQZ*TQZ+DWQZ*WQZ+ 

& DTTT*TTT+DWWW*WWW+DQQQ*QQQ+DZZZ*ZZZ 

C 

L24- DTTW*M12*M12*M22 + DTWW*M12*M22*M22 + 

& DTTQ*M12*M12*M32 + DTQQ*M12*M32*M32 + 

& DTTZ*M12*M12*M42 + DTZZ*M12*M42*M42 + 

& DWWQ*M22*M22*M32 + DWQQ*M22*M32*M32 + 

& DWWZ*M22*M22*M42 + DWZZ*M22*M42*M42 + 

& DQQZ*M32*M32*M42 + DQZZ*M32*M42*M42 + 

& DTWQ*M12*M22*M32 + DTWZ*M12*M22*M42 + 

& DTQZ*M12*M42*M32 + DWQZ*M22*M42*M32 + 

& DTTT*M12*M12*M12 + DWWW*M22*M22*M22 + 

& DQQQ*M32*M32*M32 + DZZZ*M42*M42*M42 

C 

L31=L21 
L32=L22 
L33-L23 
L34>L24 

DEFINITION OF ALFA(I) AND BETA(I) 


D1 »N32*L25 + N33*L35 
D2 »N32*L26 + N33*L36 
D3 =N32*L27 + N33*L37 
C 

Dll=-P 

D12=OMEGAO 

D21--2*OMEGAO 

D22=-P 

D23-=2*OMEGAO 

D32--OMEGAO 

D33--P 

C 

ALFA2=(D2-D21*D1/D11-D23*D3/D33)/ 

& (D22-D21*D12/D11-D23*D32/D33) 

ALFA1=(D1-D12*ALFA2)/Dll 
ALFA3-(D3-D32*ALFA2)/D33 
C 

D1 =N42*L25 + N43*L35 
D2 -N42*L26 + N43*L36 
D3 =-N42*L27 + N43*L37 
C 

D11--Q 

D12-OMEGAO 

D21--2*OMEGAO 
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D22--Q 

D23-2*OMEGAO 

D32--OMB6AO 

D33--Q 

C 

BETA2-(D2-D21*D1/D11-D23*D3/D33)/ 

& (D22-D21*D12/D11-D23*D32/D33) 

BETAl-(D1-D12*BETA2)/Dll 
BETA3-(D3-D32*BETA2)/D33 
C 

L21A-2*C11* (ALFA1*M31*M33+BETA1*M31*M34) + 
C12* (ALFAl* (M21*M33+M23*M31) + 

BETAl* (M21*M34+M24*M31) ) 

L22A-2*C11*(ALFAl*M3 2 *M3 3+BETA1*M3 2 *M3 4) 
+2*C11* (ALFA3*M31*M33+BETA3*M31*M34) 
+C12*(ALFAl*(M22*M33+M32*M23) 

+BETA1*(M22*M34+M24*M32)) 

+C12* (ALFA3* (M21*M33+M23*M31) 

+BETA3* (M21*M34+M24*M31) ) 

L23A-2*C11*( ALFA2 *M31 *M3 3 +BETA2 *M31 *M3 4) 
+2 * Cl1*(ALFA3 *M3 2 *M3 3 +BETA3 *M3 2 *M3 4) 
+C12* (ALFA2* (M21*M33+M23*M31) 

+BETA2* (M21*M34+M24*M31) ) 

+C12*(ALFA3*(M22*M33+M23*M32) 

+BETA3* (M22*M34+M24*M32) ) 

L2 4A-2 * Cl 1 * (ALFA2 *M3 2 *M3 3 +BETA2 *M3 2 *M3 4) 
+C12*(ALFA2*(M22*M33+M23*M32) 
+BETA2* (M22*M34+M24*M32) ) 

L31A=2 *C22 * (ALFAl*M31 *M33+BETA1 *M31 *M34) 
+C21*(ALFAl*(M21*M33+M23*M31) 

+BETA1* (M21*M34+M24*M31) ) 


L32A=2*C22*( ALFAl *M3 2 *M3 3 +BETA1 *M3 2 *M3 4) 
+2 * C2 2 *(ALFA3 *M31*M3 3 +BETA3 *M31*M3 4) 
+C21*(ALFAl*(M22*M33+M32*M23) 

+BETA1* (M22*M34+M24*M32) ) 

+C21*(ALFA3*(M21*M33+M23*M31) 

+BETA3* (M21*M34+M24*M31) ) 

C 
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L33A-2*C22* (ALFA2*M31*M33+BETA2*M31*M34) 
& +2*C22* (ALFA3*M32*M33+BETA3*M32*M34) 

& +C21* (ALFA2* (M21*M33+M23*M31) 

& +BETA2* (M21*M34+M24*M31) ) 

& +C21* (ALFA3* (M22*M33+M23*M32) 

& +BETA3* (M22*M34+M24*M32)) 


L3 4A-2*C22 * (ALFA2 *M32 *M3 3 +BETA2 *M3 2 *M34) 

& +C21* (ALFA2* (M22*M33+M23*M32) 

& ■^BETA2* (M22*M34+M24*M32) ) 

C 

L21=L21A+L21*B1*U*U-A13*ZGB*M11**3/6.D0 
L22-L22A+L22*B1*U*U-A13*ZGB*M12*M11**2/2.D0 
L23=L23A+L23*B1*U*U-A13*ZGB*M11*M12**2/2.D0 
L24-L24A+L24*B1*U*U-A13*ZGB*M12**3/6.D0 
L31-L31A+L31*B2*U*U-A23*ZGB*M11**3/6.D0 
L32»L32A+L32*B2*U*U-A23*ZGB*M12*M11**2/2.D0 
L33-L33A+L33*B2*U*U-A23*ZGB*M11*M12**2/2.D0 
L34-L34A+L34*B2*U*U-A23*ZGB*M12**3/6.D0 
L41--0.5*M11*M11*(M21-U*Mll/3.0) 

L42--Mil* (M12*M21+0.5*M11*M22 - 0.5*U*M11*M12) 
L43--M12*(Mll*M22+0.5*M12*M21-0.5*U*M11*M12) 
L44=-0.5*M12*M12*(M22-U*M12/3.0) 

R11=N12 *L21+N13*L31+N14*L41 
R12=N12*L22+N13*L32+N14*L42 
R13=N12*L23+N13*L33+N14*L43 
R14=N12*L24+N13*L34+N14*L44 
R21=N22*L21+N23*L31+N24*L41 
R22-N22*L22+N23*L32+N24*L42 
R23-N22*L23+N23*L33+N24*L43 
R24=N22*L24+N23*L34+N24*L44 

EVALUATE DALPHA AND DOMEGA 

UINC=0.001 
UR =U+UINC 
UL -U-UINC 
U =UR 

A(1,1)-0.0D0 
A(1,2)=0.0D0 
A(1,3)-1.0D0 
A(1,4)-0.0D0 

A{2,1)-A13*ZGB+B1*U*U*K1 
A(2,2)»A11*U +B1*U*U*K2 

A(2,3)-A12*U +B1*U*U*K3 

A(2,4). B1*U*U*K4 

A(3,1)=A23*ZGB+B2*U*U*K1 
A(3,2)=A21*U +B2*U*U*K2 
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A{3,3)-A22*U +B2*U*U*K3 

A(3,4)» B2*U*U*K4 

A(4,l)--U 
A(4,2)-1.0D0 
A(4,3)-0.0D0 
A(4,4)-0.0D0 
C 

CALL RG(4,4,A.WR,WI,0,YYY,IV1,FV1,IERR) 

CALL DSTABL(DECS,WR,WI,FREQ) 

ALPHR-DEOS 
OMEGRb FREQ 
C 

U«UL 

C 

A(1,1)-0.0D0 
A(1,2)-0.0D0 
A(1,3)-1.0D0 
A(l,4)-O.ODO 
A(2,l)«A13*ZGB+B1*U*U*K1 
A(2,2)=A11*U +B1*U*U*K2 

A(2,3)«A12*U +B1*U*U*K3 

A(2,4)* B1*U*U*K4 

A(3,1)=A23*ZGB+B2*U*U*K1 
A(3,2)=A21*U +B2*U*U*K2 

A(3,3)-A22*U +B2*U*U*K3 

A(3,4)- B2*U*U*K4 

A(4,l)=-U 
A(4,2)=1.0D0 
A(4,3)«0.0D0 
A(4,4)»0.0D0 
C 

CALL RG(4,4,A,WR,WI,0,YYY,IVl,FVl,lERR) 

CALL DSTABL(DECS,WR,WI,FREQ) 

ALPHL-DEOS 

OMEGL*=FREQ 

C 

DALPHA=(ALPHR-ALPHL)/(UR-UL) 

DOMEGA-(OMEGR-OMEGL)/(UR-UL) 

EVALUATION OF HOPF BIFURCATION COEFFICIENTS 

COEFl-3.0*Rll+R13+R22+3.0*R24 
COEF2-3.0*R21+R23-R12-3.0*R14 
AMU2 =-COEFl/(8.0*DALPHA) 

BETA2=0.25*COEF1 
C TAU2 =-(C0EF2-DOMEGA* COEFl/DALPHA)/(8.0 *OMEGAO) 

PER =2.0*3.1415927/OMEGA0 
PER -PER*U/L 

WRITE (20,2001) TC*U0/L,COEFl,DALPHA, 

& PER,AMU2,TAU2,DOMEGA 

1 CONTINUE 
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STOP 

1001 FORMAT (' ENTER NUMBER OF DATA LINES') 

1002 FORMAT (' ENTER UO, ZG, AND DSAT') 

1003 FORMAT {' ENTER BOW PLANE TO STERN PLANE RATIO') 

1004 FORMAT (' ENTER ZG') 

2001 FORMAT (7E14.5) 

2002 FORMAT (6E14.5) 

2003 FORMAT (9F11.5) 

3001 FORMAT (215) 

END 

SUBROUTINE DSOMEG(UK,WR,WI,OMEGA,CHECK) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION WR(4),WI(4) 

CHECK--1.0D+25 
DO 1 1-1,4 

IF (WR(I).LT.CHECK) GO TO 1 
CHECK-WR(I) 

IJ-I 

1 CONTINUE 

OMEGA-DABS (WI (IJ) ) 

IF (WI(IJ).GT.O.DO) IJK-IJ 
IF (WI(IJ).LT.O.DO) IJK-IJ-1 
RETURN 
END 

SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION WR(4),WI(4) 

DEOS--1.0D+20 
DO 1 1-1,4 

IF (WR(I).LT.DEOS) GO TO 1 
DEOS-WR(I) 

IJ-I 

1 CONTINUE 
OMEGA-WI (IJ) 

OMEGA-DABS(OMEGA) 

RETURN 

END 

SUBROUTINE NORMAL(T) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION T(4,4),TN0R(4,4) 

CFAC-T(1,1)**2+T(l,2)**2 

IF (DABS(CFAC).LE.(l.D-10)) STOP 4001 

TNOR(1,1)=1.DO 

TNOR(2,l)=(T(l,l)*T(2,1)+T(2,2)*T(1,2))/CFAC 
TNOR(3,l)-(T(l,l)*T(3,l)+T(3,2)♦T(l,2))/CFAC 
TNOR(4,l)-(T(l,l)*T(4,l)+T(4,2)*T(l,2))/CFAC 
TNOR(1,2)-O.DO 

TNOR(2,2)-(T(2,2)*T(l,l)-T(2,1)*T(1,2))/CFAC 
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TNOR(3,2)-(T(3,2)*T(l,l)-T(3,1)*T(1,2))/CFAC 
TN0R(4,2)-(T(4,2)*T(1,1)-T(4,1)*T(1,2))/CFAC 
DO 1 1-1,4 
DO 2 J-1,2 

T(I,J)-TNOR(I.J) 

2 CONTINUE 
1 CONTINUE 
RETURN 
END 


SUBROUTINE MULT{TINV,A,T,A2) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION TINV(4,4),A(4.4),T(4,4),A1(4,4),A2(4,4) 
DO 1 1-1,4 
DO 2 J-1,4 
A1(I,J)-O.DO 
A2 (I, J) -O.DO 

2 CONTINUE 
1 CONTINUE 

DO 3 1-1,4 
DO 4 J-1,4 
DO 5 K-1,4 

Aid, J)-A(I,K)*T{K, J)-*-Al(I, J) 

5 CONTINUE 
4 CONTINUE 

3 CONTINUE 
DO 6 1-1,4 

DO 7 J-1,4 
DO 8 K=l,4 

A2(I,J)-TINV(I,K)*A1(K,J)+A2(I,J) 

8 CONTINUE 

7 CONTINUE 

6 CONTINUE 
DO 11 1-1,4 

C WRITE (*,101) (A(I,J),J-1,4) 

11 CONTINUE 
DO 12 1=1,4 

C WRITE (*,101) (T(I,J),J-1,4) 

12 CONTINUE 
DO 10 1=1,4 

C WRITE (*,101) (A2(I,J),J-1,4) 

10 CONTINUE 

C WRITE (*,101) A2(l,l) 

RETURN 

101 FORMAT (4E15.5) 

END 
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APPENDIX C - ADAMS-BASHFORTH SIMULATION PROGRAM 

PROGRAM ABSIMl.FOR 

SUBOFF SIMULATION PROGRAM USING A FOURTH ORDER 
ADAMS-BASHFORTH INTEGRATION SCHEME 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

INTEGER LDA,NA 
PARAMETER (NA»4,LDA»4) 

DOUBLE PRECISION L,MW,MASS. IY,MQDOT,MWDOT,MQ, ICP, LEN, 1 .7, 

* MDS,MDB,MDELTA.DELTA,U0,TII Z.TC.Kl K2, 

* K3,K4,F1(4) ,F2 (4) ,F3 (4) ,F^ (4) , 

* A(LDA,NA), -SAT 
COMPLEX EVAL(NA) 

DIMENSION XL(25),BR(25),VEC1(25),VEC2(25) 

OPEN (10,FILE='SIM.DAT',STATUS*'OLD') 

OPEN (30,FILE*'ANG.DAT',STATUS*'OLD') 

OPEN (20,FILE*'SIM.RES',STATUS*'NEW') 

ENTER SPEEDS AND TIME DATA 

READ(10,*) UO,U,TC,TSIM,DELT,IPRNT,ALPHA 
U*U*U0 

GEOMETRIC PROPERTIES AND HYDRODYNAMIC COEFFICIENTS 

DSAT *0.4000D0 

PI =4.0*DATAN(1.0D0) 

RHO »1.94D0 

CDZ -0.5000D0 

L =13.979200 

WEIGHT*1556.2363D0 
lY *561.3200 

MASS =WEIGHT/32.2 
TC *TC*L/U0 

ZQDOT=-6.3300E-04*0.5*RHO*L**4 
ZWDOT=-1.4529E-02*0.5*RHO*L**3 
ZQ * 7.5450E-03*0.5*RHO*L**3 
ZW =-1.3910E-02*0.5*RHO*L**2 
ZDS =(-5.6030E-03*0.5*RHO*L**2) 

ZDB *(-5.6030E-03*0.5*RHO*L**2) 
MQDOT=-8.8000E-04*0.5*RHO*L**5 
MWDOT=-5.6100E-04*0.5*RHO*L**4 
MQ *-3.7020E-03*0.5*RHO*L**4 
MW = 1.0324E-02*0.5*RHO*L**3 
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MDS -(-2.4090E-03*0.5*RHO*L**3) 

MDB -( 2.4090E-03*0.5*RHO*L**3) 

C 

2G -0.4D0 

ZB -0.ODO 

XG -0.ODO 

XB -O.ODO 

ZGB -ZG-ZB 

ZD» ZDS +ALPHA* ZDB 
MD-MDS +ALPHA*MDB 

SET INITIAL CONDITIONS 

THETA*0.ODO 
READ(30,*) THETA 
W »0.1D0 

Q =0.0D0 

Z =0.0D0 

Z =Z*L 

DETERMINE [A] AND [B] COEFFICIENTS 

DV *(MASS-ZWDOT)*(lY-MQDOT) 

& -(MASS*XG+2QDOT)*(MASS*XG+MWDOT) 

A11DV=(lY-MQDOT)*ZW+(MASS*XG+ZQDOT)*MW 
A12DV* (lY-MQDOT) * (MASS+ZQ) + (I4ASS*XG+ZQDOT) * (MQ-MASS*XG) 
A13DV«-(MASS*XG+ZQDOT)*WEIGHT 
BIDV =(lY-MQDOT)*ZD+(MASS*XG+ZQDOT)*MD 
A21DV=(MASS-ZWDOT)*MW+{MASS*XG+MWDOT)*ZW 
A22DV*(MASS-ZWDOT)*(MQ-MASS*XG) 

& +(MASS*XG+MWDOT)♦(MASS+ZQ) 

A23DV=-(MASS-ZWDOT)*WEIGHT 
B2DV =(MASS-ZWDOT)*MD+(MASS*XG+MWDOT)*ZD 

A11=A11DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 
B1 =B1DV/DV 
B2 =B2DV/DV 

CALCULATE THE CONTROL GAINS 

POLE=1.0/TC 
ALPHA3»4.0*POLE 
ALPHA2-6.0*POLE**2 
ALPHAl-4.0*POLE**3 
ALPHAO- POLE**4 
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K4-AL.PHA0/ { (B1*A21-B2*A11) *U0**4 
& +{B1*A23-B2*A13)*ZGB*U0**2) 

A2M-B1*U0**2 

A3M-B2*U0**2 

A0M=-(A11+A22)*U0-ALPHA3 
B1M-B2*U0**2 

B2M-(B2*A12-B1*A22)*U0**3 
B3M-(B1*A21-B2*A11)*U0**3 

BOM-{A11*A22-A21*A12)*U0**2-A23*ZGB-ALPHA2-B1*U0*U0*K4 
C1M-(B2*A11-B1*A21)*U0**3 
C2M-(B1*A23-B2*A13)*ZGB*U0**2 

COM-(A13*A21-A23*A11)*ZGB*U0+ALPHA1 
& -(B2+B1*A22-B2*A12)*K4*U0**3 

K2=C1M*B0M*A3M-B1M*C0M*A3M-C1M*B3M*A0M 
K2-K2/(C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

Kl-(C0M-C2M*K2)/CIM 
K3= (AOM-A2M*K2 ) /A3M 

WRITE(*,*) K1,K2,K3,K4 

DETERMINE THE EIGENVALUES 

A{1,1)=0.0D0 

A(1,2)-0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2,l)=A13*ZGB+B1*U*U*K1 
A(2,2)-A11*U +B1*U*U*K2 

A(2,3)=A12*U +B1*U*U*K3 

A(2,4)= B1*U*U*K4 

A(3,1)-A23*ZGB+B2*U*U*K1 
A(3,2)=A21*U +B2*U*U*K2 

A(3,3)=A22*U +B2*U*U*K3 

A{3,4)= B2*U*U*K4 

A(4,l)--U 
A(4,2)=1.0D0 
A(4,3)«0.0D0 
A(4,4)=0.0D0 

CALL DEVLRG(NA,A,LDA,EVAL) 

CALL DWRCRN ('EVAL',1,NA,EVAL,1,0) 

DEFINE THE LENGTH X AND BREADTH B TERMS FOR THE 
INTEGRATION 

XL( 1)= O.OOOODO 
XL{ 2)= O.IOOODO 
XL( 3)- 0.2000D0 

XL( 4)= 0.3000D0 
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XL( 5)- 0.4000D0 

XL( 6)- 0.5000D0 

XL( 7)- 0.6000D0 

XL( 8)- 0.7000D0 

XL( 9)- l.OOOODO 
XL(10)- 2.0000D0 

XL (ID- 3.0000D0 
XL(12)- 4.0000D0 

XL(13)- 7.7143D0 

XL(14)- lO.OOOODO 
XL(15)- 15.1429D0 
XL(16)- 16.0000D0 
XL(17)- 17.0000D0 
XL(18)- 18.0000D0 
XL(19)- 19.0000D0 
XL(20)- 20.0000D0 
XL(21)= 20.1000D0 
XL(22)= 20.2000D0 
XL(23)= 20.3000D0 
XL(24)= 20.4000D0 
XL(25)= 20.4167D0 

DO 102 N = 1,25 

XL(N) = (L/20)*XL(N) - L/2 
102 CONTINUE 


BR( 1)= O.OOODO 
BR( 2)= 0.485D0 
BR( 3)- 0.658D0 
BR( 4)= 0.778D0 
BR( 5)= 0.871D0 
BR( 6)= 0.945D0 
BR( 7)= I.OIODO 
BR( 8)= 1.060D0 
BR( 9)= 1.180D0 
BR(10)= 1.410D0 
BR(11)= 1.570D0 
BR(12)= 1.660D0 
BR(13)- 1.670D0 
BR(14)= 1.670D0 
BR(15)= 1.670D0 
BR(16)= 1.630D0 
BR(17)= 1.370D0 
BR(18)- 0.919D0 
BR(19)- 0.448D0 
BR(20)- 0.195D0 
BR(21)- 0.188D0 
BR(22)- 0.168D0 
BR(23)= 0.132D0 
BR(24)= 0.053D0 
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BR(25)= O.OOODO 
C 

PISIM -TSIM/DELT 
ISIM -PISIM 
C 

C SIMULATION BEGINS 

C 

DO 100 I-1,ISIM 
C 

C CALCULATE THE DRAG FORCE, INTEGRATE THE DRAG OVER THE 

C VEHICLE 

C 

DO 101 K-1.25 
UCF-W-XL(K)*Q 

VECl(K)=BR(K)*UCF*DABS(UCF) 

VEC2(K)-BR(K)*UCF*DABS(UCF)*XL(K) 

101 CONTINUE 

CALL TRAP(25,VECl,XL,HEAVE) 

CALL TRAP(25,VEC2,XL,PITCH) 

HEAVE-0.5*RHO*CDZ*HEAVE 
PITCH-0.5*RHO*CDZ*PITCH 
C 

C COUPLING EQUATIONS 

C 

DWDV-(lY-MQDOT)*HEAVE+(MASS*XG+ZQDOT)*PITCH 
DQDV-(MASS-ZWDOT)*PITCH+(MASS*XG+MWDOT)*HEAVE 
CIDV-(lY-MQDOT)*MASS*ZG*Q**2 
& -(MASS*XG+ZQDOT)*MASS*ZG*W*Q 

C2DV--(MASS-ZWDOT)*MASS*ZG*W*Q 
& +(MASS*XG+MWDOT)*MASS*ZG*Q**2 

C 

DW-DWDV/DV 

DQ-DQDV/DV 

Cl-ClDV/DV 

C2-C2DV/DV 

C 

C LIMITING CONDITION 

C 

IF (ABS(THETA) .GT. PI/4) THEN 

WRITE (*,*) 'BUBBLE EXCEEDED +/- 45, BOAT WAS LOST' 
STOP 
ENDIF 
C 

C CONTROL EQUATIONS 

C 

DELTA = (K1*THETA+K2*W+K3*Q+K4*Z) 

DELTA = DSAT*DTANH(DELTA/DSAT) 

C IF (DELTA .GT. 0.4) DELTA-0.4 

C IF (DELTA .LT. -0.4) DELTA--0.4 

C 

THEDOT-Q 
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WDOT -A11*U*W+A12*U*Q+A13*ZGB*DSIN(THETA) 
& +B1*U*U*DELTA+DW+C1 

QDOT -A21*U*W+A22*U*Q+A23*ZGB*DSIN(THETA) 
& +B2*U*U*DELTA+DQ-*-C2 

ZDOT --U*DSIN(THETA)+W*DCOS(THETA) 

IF (I.GT.3) GO TO 150 

INITIAL FIRST ORDER INTEGRATION 

THETA - THETA + DELT*THEDOT 

W - W + DELT*WDOT 

Q - Q + DELT*QDOT 

Z - Z + DELT*ZDOT 

Fl(I)»THEDOT 
F2(I)-WDOT 
F3(I)-QDOT 
F4(I)-ZDOT 

ADAMS-BASHFORTH INTEGRATION 
150 F1(4)=Q 

F2(4)-A11*U*W+A12*U*Q+A13*ZGB*DSIN(THETA) 

& +B1*U*U*DELTA+DW+C1 

F3(4)=A21*U*W+A22*U*Q+A23*ZGB*DSIN(THETA) 

& +B2*U*U*DELTA+DQ+C2 

F4(4)=-U*DSIN(THETA)+W*DCOS(THETA) 

THETA-THETA+(DELT/24.0)*(55.0*F1(4) 

1 -59.0*F1(3)+37.0*F1(2)-9.0*F1(1) ) 

W =W+(DELT/24.0)*(55-0*F2(4) 

1 -59.0*F2(3)+37.0*F2(2)-9.0*F2(1)) 

Q -Q+(DELT/24.0)*(55.0*F3(4) 

1 -59.0*F3(3)+37.0*F3(2)-9.0*F3(1)) 

Z =Z+(DELT/24.0)*(55.0*F4(4) 

1 -59.0*F4(3)+37.0*F4(2)-9.0*F4(1) ) 

F1(1)-F1(2) 

F1(2)=F1(3) 

F1(3)-F1(4) 

F2(1)=F2(2) 

F2(2)-F2(3) 

F2(3)=F2(4) 

F3(1)-F3(2) 

F3(2)»F3(3) 

F3(3)-F3(4) 

F4(l)-F4(2) 

F4(2)-F4(3) 

F4(3)-F4(4) 

TIME-I*DELT 
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J-J+1 

IF (J.NE.IPRNT) GO TO 100 
ALFA-W/U 

ALFA«(-DATAN(ALFA))*180/PI 

WRITE (20,991) TIME*U0/L,Z/L,THETA*180/PI, 

& DELTA,W,Q,ALFA 

WRITE (*,991) TIME*U0/L,Z/L,THETA*180/PI, 

& DELTA,W,Q,ALFA 

J-0 
C 

100 CONTINUE 
STOP 

991 FORMAT (2X,F8.2,6(2X,F8.4)) 

END 

C 

SUBROUTINE TRAP(N,A,B,OUT) 

C 

C NUMERICAL INTEGRATION ROUTINE USING THE TRAPEZOIDAL RULE 
C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION A(l),B(1) 

N1=N-1 
OUT=0.0 
DO 1 1=1,N1 

OUT1=0.5* (A (I) +A(I + 1) ) * (B(I + 1) -Bd) ) 

OUT =OUT+OUTl 
1 CONTINUE 
RETURN 
END 
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APPENDIX D - ROOT LOCUS PROGRAM 

% THIS MATLAB PROGRAM FINDS THE ROOT LOCUS FOR A GIVEN UO, Zg, 
% Alpha, AND Tc. 

W-1556.2363;Bl-1556.2363; 

Iy-561.32; 
g«32.2; 
m-W/g; 
rho-1.94; 

L-13.9792; 
xg*0;zg».4; 
zgb=.4; 

ndl=.5*rho*L^2; 
nd2=.5*rho*L*3; 
nd3=. 5*rho*L'^4 ; 
nd4=.5*rho*L^5; 

e * zeros(100,4); 

X » zeros(4,1); 

Y * zeros(4,1); 

Z = zeros(4,1); 

Alpha=0; 

Zqdnd=-6.33e-4;Zwdnd=-1.4529e-2;Zqnd*7.545e-3;Zwnd=-1.391e-2; 
Zds=.5*(-5,603e-3);Zdb=-.5*5.603e-3;Zdltnd=(Zds+Alpha*Zdb); 
Mqdnd=-8.8e-4;Mwdnd=-5.61e-4;Mqnd=-3.702e-3;Mwnd=l.0324e-2; 
Mds=.5*(-0.002409);Mdb=.5*0.002409;Mdltnd=(Mds+Alpha*Mdb); 

Zqd=nd3 * Zqdnd;Zwd=nd2 * Zwdnd;Zq=nd2 * Zqnd;Zw=ndl* Zvmd; 

Zdlt=ndl* Zdltnd; 

Mqd=nd4 *Mqdnd;Mwd=nd3 *Mwdnd;Mq=nd3 *Mqnd;Mw=nd2 *Mwnd; 

Mdlt=nd2 *Mdltnd; 

Dve(m-Zwd)*(ly-Mqd)-(m*xg+Zqd)*(m*xg+Mwd); 

allDv=(ly-Mqd)*Zw+(m*xg+Zqd)*Mw; 

al2Dv= (ly-Mqd) * (m+Zq) + (m*xg+Zqd) * (Mq-in*xg) ; 

al3Dv=-(m*xg+Zqd)*W; 

blDv=(ly-Mqd)*Zdlt+(m*xg+Zqd)*Mdlt; 

a21Dv=(m-Zwd)*Mw+(m*xg+Mwd)*Zw; 

a22Dv» (m-Zwd) * (Mq-m*xg) + (m*xg+Mwd) * (m+Zq) ; 

a23Dv=-(m-Zwd)*W; 

b2Dv»(m-Zwd)*Mdlt+(m*xg+Mwd)*Zdlt; 

all=allDv/Dv;al2«al2Dv/Dv;al3=al3Dv/Dv; 
a21=a2IDv/Dv;a2 2 »a22Dv/Dv;a2 3=a2 3Dv/Dv; 
bl=blDv/Dv;b2=b2Dv/Dv; 
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u - 9; 

for j ■ 1:240 
U - j/20; 
sl(j)- U; 

Fn-sqrt (U‘^2/ (zgb*g)) ; 
tconst-4.7512; 
pole--l/{tconst*L/u); 

pl»[pole pole-0.01 pole-0.02 pole-0.03]; 
a«[0 0 1 0;al3*zgb all*u al2*u 0;a23*zgb a21*u a22*u 0;... 
-u 1 0 0] ; 

aU=[0 0 1 0;al3*zgb all*U al2*U 0;a23*zgb a21*U a22*U 0;-U 
10 0 ]; 


b«[0;bl*u*2;b2*u^2;0]; 
bU-[0;bl*U^2;b2*U^2;0]; 

Kl=place(a,b,pi); 

Y = eig(aU-bU*Kl); 
e(j, :) = Y' ; 

end 

for k = 1:200 

U = 8.9 + k/1000; 
s2(k) = U; 

Fn=xsqrt (U'^2/ (zgb*g)) ; 
tconst=4.7512; 
pole=-l/(tconst*L/u); 

pl=[pole pole-0.01 pole-0.02 pole-0.03]; 
a=[0 0 1 0;al3*zgb all*u al2*u 0;a23*zgb a21*u a22*u 0;... 
-u 1 0 0]; 

aU=[0 0 1 0;al3*zgb all*U al2*U 0;... 
a23*zgb a21*U a22*U 0;-U 100]; 

b= [0;bl*u*2;b2*u^2;0]; 
bU=[0;bl*U^2;b2*U^2;0]; 

Kl=place{a,b,pl); 

Z = eig(aU-bU*Kl); 
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f (k, ;) - Z' ; 
end 


cucis ([- .4 .1 - .8 .8]) 
plot(real(e(1,1)),imag(e(1,1)),'x',... 
real(e{l,2)),imag(e(l,2)),'x',... 
real (e (1,3)), imag (e (1,3)), 'x',... 
real(e(1,4)),imag(e(1,4)),'x') 

hold 

plot(real(e(:,1)),imag(e(:,1). 
real(e(:,2)),imag(e(:,2). 
real(e(:,3)),imag(e(:,3). 
real(e(:,4)),imag(e(:,4)),'.') 
plot(real(f(:,1)),imag(f(:,1)),'.',... 
real(f(:,2)),imag(f(:,2)),'.',... 
real(f(:,3)),imag(f(:,3)),'.',... 
real(f(:,4)),imag(f(:,4)),'.') 
grid,title('EIGENVALUES FOR SPEEDS 1 TO 100 FPS') 
%title('DETAILED LOW SPEED ROOT LOCUS PLOT') 
pause 
hold 
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